

%% load everything of interest
load('data.mat')
load('sim_eq_round1.mat')
load('parameter_estimates_round1_replication.mat','specif')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% normalize continuous variables
% data.gdp_pc = (data.gdp_pc-mean(data.gdp_pc))/std(data.gdp_pc);
% data.pop = (data.pop-mean(data.pop))/std(data.pop);
% data.polity = (data.polity-mean(data.polity))/std(data.polity);
data.terrain = (data.terrain-mean(data.terrain))/std(data.terrain);
dd = [data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
dd = (dd-mean(dd(:)))/std(dd(:));
data.USdist = dd(:,1);
data.UKdist = dd(:,2);
data.FRdist = dd(:,3);
data.RUdist = dd(:,4);
data.CHdist = dd(:,5);
clear dd


%% estimates of mean utilities
SV=3000; % # of starting values
X=[ones(N,1) data.terrain]; K=size(X,2);
Z_US=data.USdist; KZ=size(Z_US,2);
Z_UK=data.UKdist;
Z_France=data.FRdist;
Z_Russia=data.RUdist;
Z_China=data.CHdist;
Z=[Z_US,Z_UK,Z_France,Z_Russia,Z_China];
reb={1:K,[]}; mpow={1,1:KZ};
Y=arrayfun(@(n)find(sum(YY==table2array(outcomes(n,4:end)),2)==size(YY,2)),1:N)';

options=optimset('Display','off','HessFcn',@(x,lambda)HessPhat(x,lambda,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));

load('parameter_estimates_round1_replication.mat','srngSV')
rng(srngSV)
ltheta=length(reb{1})+4+length(mpow{1})+5+length(reb{2})+length(mpow{2})+10+size(DMES{1,1},2);
theta0=[zeros(ltheta,1),randn(ltheta,SV-1)/10]; % starting values
theta=theta0; 
l0=1:SV; % maximized likelihood values across starting values
exitn=1:SV; % Knitro exit flags across starting values
f0 = exitn;

time=0;
parfor sv=1:SV
    tic
    f0(sv) = Phat(theta0(:,sv),Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
    if ~isfinite(f0(sv))
        theta(:,sv) = NaN(length(theta0(:,sv)),1);
        exitn(sv) = NaN;
        l0(sv) = NaN; 
        time=time+toc;
    else
        [theta(:,sv),l0(sv),exitn(sv)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
            theta0(:,sv),[],[],[],[],[],[],[],[],options,'knitro1.opt'); % use knitro0.opt to check derivatives
        time=time+toc;
    end
end


%%
clearvars -except exitn K l0 N SV theta theta0 time X Y YY Z reb mpow srngSV specif
save('parameter_estimates_round1.mat')


%% verify optimality
[~,use] = min(l0);
thetaStar = theta(:,use);
load('parameter_estimates_round1_replication.mat')
[~,use] = min(l0);
thetaStar0 = theta(:,use);
disp(' ')
disp('verify round 1 results:')
disp(['max. rel. diff. between estimates = ',num2str(max(abs(thetaStar0-thetaStar)./abs(thetaStar0)))])



